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ABSTRACT 


Multiple  Model  Adaptive  Estimation  (MMAE)  is  a  Bayesian  technique  that  applies 
a  bank  of  Kalman  filters  to  predict  future  observations.  Each  Kalman  filter  is  based  on  a 
different  set  of  parameters  and  hence  produces  different  residuals.  The  likelihood  of  each 
Kalman  filter’s  prediction  is  determined  by  a  magnitude  of  the  residuals.  Since  some 
researchers  have  obtained  good  forecasts  using  a  single  Kalman  filter,  we  tested 
MMAE’s  ability  to  make  time  series  predictions.  Our  Kalman  filters  have  a  dynamics 
model  based  on  a  Box-Jenkins  Auto-Regressive  Moving  Average  (ARMA)  model  and  a 
measure  model  with  additive  noise.  The  time-series  prediction  is  based  on  the 
probabilistic  weighted  Kalman  filter  predictions.  We  make  a  probability  interval  about 
that  estimate  also  based  on  the  filter  probabilities.  In  a  Monte  Carlo  analysis,  we  test  this 
MMAE  approach  and  report  the  results  based  on  many  different  criteria.  Our  analysis 
tests  the  robustness  of  the  approach  by  testing  its  ability  to  make  predictions  when  the 
Kalman  filter  dynamics  models  did  not  match  the  data  generation  time-series  model.  Our 
analysis  indicates  benefits  in  applying  multiple  model  adaptive  estimation  for  time  series 
analysis. 


MULTIPLE  MODEL  ADAPTIVE  ESTIMATION 


FOR  TIME  SERIES  ANALYSIS 

CHAPTER  I:  INTRODUCTION 


1.1  Overview 

In  this  chapter  we  introduce  the  development  of  Multiple  Model  Adaptive 
Estimation  (MMAE)  [4, 5, 11, 13, 14, 16,  21, 22]  for  Time  Series  Analysis.  In  Section 
1.2,  we  present  the  problem  perspective  and  background  for  time  series  analysis.  In 
Sections  1.3  and  1.4,  we  define  the  details  of  research  objectives,  motivation  for  using 
MMAE  for  time  series  analysis,  scope  of  research,  limitations  and  assumptions  to 
accomplish  this  research.  Section  1.5  presents  the  methodology.  Finally  Section  1.6  is  the 
overview  of  following  chapters. 

1.2  Problem  Perspective 

A  time  series  is  set  of  observations  on  a  quantitative  variable  collected 
sequentially  in  time  [1,  3,  15].  Many  sets  of  data  appear  as  time  series:  a  monthly  amount 
of  goods  shipped  from  a  factory,  monthly  aircraft  spare  part  stock  levels  in  an  air  force, 
weekly  numbers  of  accidents  on  a  road,  the  daily  closing  values  of  a  stock  market,  hourly 
yields  of  a  chemical  reaction  and  many  other  examples.  We  can  find  examples  of  time 
series  from  many  different  fields:  economics,  engineering,  natural  or  social  sciences, 
military,  business  and  so  on. 


In  business,  people  are  always  interested  in  forecasting  future  values  of  a  time 
series  variable.  They  try  to  forecast  the  costs,  sales,  profits,  inventory,  backorders  and  etc. 
The  military  attempts  to  predict  spare  part  inventories,  aircraft  mishaps,  personnel 
requirements  and  many  other  things.  To  make  accurate  time  series  forecasts,  we  need 
time  series  analysis  techniques. 

Time  series  analysis  is  a  method  used  for  making  predictions  about  future  values 
of  a  time  series  variable  based  on  past  observations  of  that  series  and  mathematical  or 
statistical  model. 

Unlike  the  analysis  of  random  samples  of  observations  that  are  discussed  in  the 
context  of  most  other  statistics,  the  analysis  of  time  series  is  based  on  the  assumption  that 
successive  values  in  the  data  file  represent  consecutive,  dependent  measurements  taken  at 
equally  spaced  time  intervals.  The  nature  of  this  dependence  of  a  time  series  is  a  major 
interest.  Time  series  analysis  concerns  techniques  for  the  analysis  of  this  dependency. 
There  are  two  main  goals  of  time  series  analysis:  (a)  identify  the  pattern  represented  by 
the  sequence  of  observations,  and  (b)  forecast  future  values  of  the  time  series  variable. 
Once  the  pattern  is  established,  we  can  extrapolate  the  identified  pattern  to  predict  future 
events.  The  basic  structure  of  time  series  models  is: 

DATA  =  PATTERN  +  NOISE  [9] 

PATTERN:  Fundamental  relationship  describing  the  expected  or  average 

performance  of  the  data  over  time.  Implied  “assumption  of 
continuity”  of  pattern  over  forecast  period. 

NOISE:  Randomness,  "error,"  or  variability  around  the  pattern. 
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With  this  structure,  the  process  of  specifying  a  model  to  represent  a  time  series  often 
becomes  one  of  pattern  recognition. 

1.3  Research  Objectives 

The  time  series  analysis  literature  presents  many  analysis  techniques:  classical 
techniques  such  as  moving  averaging  and  smoothing  methods,  and  Box-Jenkins  ARMA, 
Autoregressive  Moving  Average,  models  are  some  of  the  most  popular  [1,  3, 9, 12,  15]. 
Classical  techniques  are  generally  appropriate  for  forecasting  time  series  that  exhibit 
stable  patterns,  which  can  be  modeled  as  relatively  straightforward  functions  of  time. 
However,  in  real-life  research  and  practice,  patterns  of  the  data  are  unclear,  individual 
observations  involve  considerable  error.  Furthermore  we  still  need  to  uncover  the  hidden 
patterns  in  the  data  and  make  good  forecasts.  The  ARMA  methodology  developed  by 
Box-Jenkins  allows  us  to  do  just  that. 

We  propose  completely  a  new  and  different  approach,  MMAE,  for  time  series 
estimation.  We  use  ARMA  models  as  our  required  dynamics  model  in  MMAE  because  of 
its  capability  of  representing  almost  all  kinds  of  time  series.  Normally,  the  application  of 
Box-Jenkins  ARMA  models  requires  to  choose  only  one  model  which  best  fits  the  data. 

In  our  new  method,  MMAE,  in  order  to  forecast  future  values  of  a  time  series,  instead  of 
choosing  one  best  ARMA  model,  we  use  all  models  -  each  model  has  different  set  of 
parameters-  that  represent  time  series  adequately.  We  basically  use  the  all  most  likely 
region  of  ARMA  parameter  space  instead  of  limiting  ourselves  to  only  one  set  of 
parameters.  MMAE  allows  us  to  run  multiple  models  and  the  Kalman  filter  structure  in 
MMAE  weights  each  model,  where  less  likely  models  receive  less  Bayesian  probability. 
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The  probabilistic  weighted  average  of  each  Kalman  filter’s  prediction  gives  the  MMAE 
forecast.  A  significant  advantage  of  our  proposed  approached  is  that  it  combines  the  Box- 
Jenkins  steps  of  identification,  and  estimating,  and  automates  it.  Our  primary  objective  is 
to  test  MMAE’s  ability  to  make  time  series  predictions. 

1.4  Scope  of  Research 

Our  research  is  basically  limited  to  the  application  of  MMAE  to  time  series 
analysis. 

MMAE  is  a  Bayesian  [9, 20,  21, 22]  technique  that  applies  bank  of  Kalman  filters 
to  predict  future  observations.  In  fact  MMAE  is  nothing  more  than  an  algorithm  that  uses 
multiple  filters  in  parallel  to  represent  the  dynamic  nature  of  the  system. 

An  overview  of  the  MMAE  algorithm  follows  [16]:  The  set-up  for  employing 
MMAE  is  to  discretize  the  continuous  space  for  each  parameter  into  a  set  of 
representative  points.  The  MMAE  algorithm  processes  measurements  (time  series  data 
values  in  this  application)  through  a  Kalman  filter  at  each  combination  of  discrete 
parameters.  Each  filter’s  residuals  determine  the  probability  of  that  filter's  parameters 
being  correct,  conditioned  on  the  measurements  processed  to  that  time.  After  processing 
all  the  available  measurements,  the  filter  probabilities  indicate  the  likelihood  of  the 
parameters  in  that  filter  being  correct  conditioned  on  the  measurements. 

The  Kalman  filter  is  simply  an  optimal  recursive,  two-step  data  processing 
algorithm  [6, 10, 12,  16,  23].  It  propagates  an  estimate  at  the  first  step  and  then  updates 
the  estimate  with  an  imprecise  measurement,  repeatedly.  Each  Kalman  filter  uses  its  own 
model  to  develop  estimation  of  current  state  independent  of  the  other  filters.  Unlike  other 
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processing  concepts,  Kalman  filter  does  not  require  all  previous  data  to  be  kept  in  storage 
and  reproduced  every  time  a  new  measurement  is  taken.  The  Kalman  filter  is  discussed  in 
detail  in  Chapter  II. 

When  we  apply  Box-Jenkins’  method,  we  should  estimate  the  parameters  of  the 
selected  ARMA  model,  by  some  sort  of  estimation  technique,  like  least  square  estimation 
or  maximum  likelihood  estimation.  However  we  may  lose  some  important  properties  of 
another  models  that  also  represent  the  time  series  adequately,  by  not  selecting  them.  We 
contend  that  if  more  than  one  model  represents  a  time  series  adequately,  they  should  each 
be  retained.  MMAE  probabilistically  weights  each  model,  incorporates  their  predictions 
based  on  their  associated  probabilities  and  makes  a  forecast.  It  has  two  major  advantages 
in  this  application  to  time  series  analysis.  First,  it  combines  the  identification  and 
estimation  steps  of  traditional  Box-Jenkins’  method.  The  entire  range  of  model 
parameters  and  model  may  be  implemented,  rather  than  selecting  a  single  model  in  the 
Box-Jenkins’  identification  stage.  Based  on  the  data,  the  algorithm  assigns  probabilistic 
weights  to  each  model.  If  two  or  more  models  appear  to  fit  the  data,  MMAE  uses  all 
model  estimates  with  appropriate  probabilistic  weights.  Second,  MMAE  adjusts  the 
model  probabilistic  weights  for  changes  in  the  data  pattern. 

1.5  Methodology 

We  use  Monte  Carlo  simulation  to  generate  notional  time  series  data,  using  an 
ARMA  model,  apply  MMAE  approach  to  this  data,  test  this  approach  and  report  the 
results  depending  on  many  different  criteria.  The  details  and  algorithm  are  given  in 
Chapter  III. 
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1.6  Thesis  Overview 


Many  researchers  have  used  MMAE  for  many  different  kinds  of  problems 
successfully  [6,  8, 11, 13, 19, 21, 24].  Furthermore  some  researchers  have  obtained  good 
forecasts  using  a  single  Kalman  filter  [10],  so  we  test  MMAE’s  ability  to  make  time 
series  predictions.  Since  we  are  unable  to  find  any  documented  application  of  MMAE  to 
time  series  analysis  in  the  literature,  we  believe  that  this  is  a  new  and  untried  approach. 

In  the  following  chapters  we  develop  this  research  in  more  detail.  Chapter  II  gives 
the  literature  review  of  relevant  subjects:  Kalman  filters,  and  Multiple  Model  Adaptive 
Estimation.  Chapter  III  presents  our  research  methodology  by  discussing  the 
development  of  MMAE  for  time  series  analysis.  Chapter  IV  gives  the  results  from  our 
analysis  based  on  many  different  criteria.  Finally  Chapter  V  summarizes  these  results  and 
concludes  with  our  recommendations  for  future  studies. 


6 


CHAPTER  II:  HISTORY  AND  BACKGROUND 


2.1  Overview 

In  this  chapter  we  present  the  background  of  Multiple  Model  Adaptive 
Estimation,  and  the  Kalman  filters.  In  section  2.2  we  discussed  history  of  MMAE  and 
Kalman  filters,  then  we  continue  with  the  overview  of  MMAE  in  Section  2.3,  Kalman 
filter  in  Section  2.4,  and  the  detailed  MMAE  algorithm  in  Section  2.5. 

2.2  History 

How  to  best  make  inferences  from  a  time  series  realization  of  data  is  a  problem 
with  a  long  history.  Much  of  the  early  work  in  time  series,  such  as  by  Bernoulli,  Gauss, 
and  Legendre,  grew  out  of  problems  in  astronomy.  The  problem  they  addressed  involved 
making  inference  as  to  the  location  of  a  “heavenly  body”  from  a  sequence  of  imperfect 
observations.  Gauss  and  Legendre  developed  and  applied  least  square  estimation  in  1795. 
R.  A.  Fisher  introduced  maximum  likelihood  estimation  in  1912.  [23] 

In  the  1800’s,  the  Danish  astronomer  T.  N.  Thiele  developed  a  recursive 
procedure  resembling  what  is  now  referred  to  as  the  Kalman  Filter  for  the  problem  of 
determining  the  distance  from  Copenhagen  to  Lund.  Kolmogorov  in  1941  and  Wiener  in 
1942  independently  developed  a  linear  minimum  mean-square  estimation  technique  that 
received  considerable  attention  and  provided  the  foundation  for  subsequent  development 
of  Kalman  Filter  theory.  [22,  23] 

Kalman’s  and  Gauss’  work  are  significantly  related.  The  Kalman  filter  can  be 
rightfully  regarded  as  an  efficient  computational  solution  of  the  weighted  least  square 
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method.  The  noise  and  initial  state  (presented  in  section  2.4)  are  essentially  assumed  by 
Kalman  to  be  Gaussian  (normally)  distributed.  Kalman  assumes  that  the  noise  is 
independent  from  one  sampling  time  to  the  next  time;  therefore  he  agrees  with  Gauss’ 
assumption  [23].  The  basic  assumptions  of  Gauss  and  Kalman  are  identical  except  that 
the  latter  allows  the  state  to  change  from  one  time  to  the  next. 

Kalman  published  his  first  famous  paper  in  1960,  and  Kalman  and  Bucy 
published  another  paper  in  1961  [23].  Almost  immediately  after  Kalman’s  paper 
appeared  (1960),  a  great  interest  from  the  aerospace  community  in  application  of 
procedures  was  aroused.  Kalman  met  with  the  individuals  at  the  NASA  Ames  Research 
Center  in  the  fall  of  1960  to  discuss  the  potential  application  of  his  work  in  the  navigation 
and  guidance  problem  for  a  possible  lunar  mission.  In  August  1961,  NASA  extended  the 
first  major  Apollo  contract  to  the  M.I.T.  Instrument  Laboratory  to  develop  an  on-board 
navigation  and  guidance  system,  which  was  to  rely  on  state-space  and  Kalman  filtering 
ideas  [22]. 

State-space  modeling  and  Kalman  filtering  is  now  the  approach  of  choice  for  a 
wide  range  of  time  series  problems  [10,  12].  Their  flexibility  in  handling  multivariate 
data  and  nonstationary  processes  provides  a  significant  advantage  over  traditional  time 
series  techniques  for  many  applications,  such  as  geophysical  exploration,  biomedicine, 
demography,  nuclear  power  plant  failure  detection,  and  macroeconomics  forecasting,  in 
addition  to  traditional  aerospace  applications. 

D.  T.  Magill  [14]  presented  the  MMAE  for  the  first  time  in  1965.  He  arranged  a 
number  of  Kalman  Filters,  each  with  different  time  invariant  plant  models  and  based  on  a 
hypothesized  parameter  realization,  and  used  the  residuals  from  these  filters  to  form  an 
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appropriately  weighted  sum  of  Kalman  filter  estimates.  He  showed  that  this  adaptive 
estimation  algorithm  produced  the  optimal  estimate  in  the  minimum  mean  square  error 
sense  for  a  Gauss-Markov  process.  Although  it  was  not  specifically  named  MMAE,  the 
structure  of  using  multiple  Kalman  Filters  was  put  into  place. 

Chang  and  Athans  extended  Magill’s  work  to  handle  discrete  systems  [4,  5]. 

Many  researchers  have  focused  on  exploiting  the  capability  of  MMAE  to  provide  state  or 
parameter  estimation  separately,  as  well  as  attempting  to  blend  them. 

Due  to  popularity  of  Kalman  filters  in  the  aerospace  industry,  USAF  frequently 
used  MMAE  to  solve  specific  Air  Force  problems.  Dr.  Peter  Maybeck  at  Air  Force 
Institute  of  Technology,  AFIT,  has  supervised  extensive  research  in  the  application  of 
MMAE  [8, 11, 13, 19,21,24], 

2.3.  Multiple  Model  Adaptive  Estimation  Review 

In  this  section  we  present  the  overview  of  MMAE.  Because  the  equations  behind 
the  MMAE  are  Kalman  filter  equations,  the  detailed  algorithm  is  presented  in  Section  2.5 
following  the  Kalman  filter  discussion. 

MMAE  is  a  Bayesian  technique  that  applies  a  bank  of  parallel  Kalman  filters  to 
predict  future  observations.  Each  Kalman  filter  is  based  on  a  different  set  of  parameters 
and  hence  produces  different  residuals.  The  likelihood  associated  with  each  Kalman  filter 
prediction  is  determined  by  the  magnitude  of  residuals. 

The  basic  structure  of  MMAE  is  shown  in  Figure  1.  The  major  strength  of 
MMAE  is  its  ability  to  reconfigure  rapidly  in  the  presence  of  failures  and  thus  provide 
accurate  state  estimates. 
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Figure  1.  Multiple  Model  Adaptive  Estimation  Algorithm  [17] 


An  overview  of  the  MMAE  algorithm  follows:  First  the  continuous  space  is 
discretized  for  each  parameter  yielding  a  set  of  representative  points.  The  MMAE 
algorithm  processes  the  measurements  through  a  Kalman  filter  at  each  combination  of 
these  discrete  parameters.  Each  filter’s  residuals  determine  the  probability  of  that  filter's 
parameters  being  correct,  conditioned  on  the  measurements  processed  to  that  time.  After 
processing  all  the  available  measurements,  the  filter  probabilities  indicate  the  likelihood 
of  the  parameters  in  that  filter  being  correct  conditioned  on  the  measurements. 

This  procedure  is  based  on  the  following  assumptions  [24]: 


to 


•  The  sampled-data  system  is  adequately  represented  by  a  linear  stochastic  state 
differential  model  for  a  given  parameter  vector  value,  resulting  in  Gaussian 
probability  density  functions,  and  can  be  described  equivalently  by  linear 
stochastic  differential  equations  [16, 17,  18].  In  case  nonlinear  models  are 
required  to  describe  the  system  adequately,  extended  Kalman  filters  [6, 14, 22] 
replace  the  linear  Kalman  filter  in  the  MMAE. 

•  The  uncertain  parameters  to  be  estimated  affect  the  system  matrices  or  the 
statistics  of  the  noises  entering  the  system. 

•  MMAE  theory  assumes  a  discrete-valued  parameter  vector.  Actually,  a  parameter 
value  typically  varies  over  a  continuous  range  of  the  parameter  space.  Thus, 
parameter  values  have  to  be  discretized  to  some  level  of  resolution  for  feasible 
implementation.  Clearly,  poor  choices  in  discrete  values  for  a  continuous 
parameter  may  result  in  poor  modeling  by  the  MMAE  elemental  filters  and  thus 
poor  estimation  from  the  entire  MMAE  itself.  This  results  when  an  elemental 
filter  within  the  MMAE  bank  does  not  have  a  good  model  of  the  system’s  current 
behaviour. 

2.4  Kalman  Filter  Theory 

In  this  section,  we  present  the  general  Kalman  filter  theory  applicable  to 
continuous  and  time-varying  systems.  The  time  series,  forecasting  application  will  be 
limited  to  discrete,  time-invariant  Kalman  filters  in  Chapter  III.  The  notation  is  based  on 
Maybeck’s  notations  [16]. 
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2.4.1  Introduction 


A  Kalman  filter  is  an  optimal,  recursive,  two-step  data  processing  algorithm.  It  is 
optimal  for  any  criterion  that  makes  sense  and  is  chosen  to  evaluate  the  performance  of  a 
system.  The  Kalman  filter  incorporates  all  available  information,  including  all 
measurements  regardless  of  their  accuracy  to  estimate  the  state  of  a  system.  The  state  of 
the  system  refers  to  a  set  of  variables  that  describe  the  inherent  properties  of  the  system  at 
a  specific  instant  in  time.  It  is  recursive,  because  the  same  two  steps  of  propagation  and 
update  are  repeated  for  each  additional  observation.  Furthermore,  it  doesn’t  require  all 
previous  data  be  kept  in  storage  and  reprocessed  every  time  a  new  measurement  is  taken. 
Filter  or  data  processing  algorithm  refers  to  a  computer  program  implements  the 
mathematical  algorithm. 

The  need  for  a  filter  is  apparent  from  the  Figure  2.  Often  times  the  state  of  system 
can’t  be  measured  directly,  and  measuring  devices  are  used  to  infer  these  values. 
Unfortunately  any  measurement  contains  some  degree  of  noise,  bias  and  device 
inaccuracies,  and  a  valuable  estimate  must  be  derived  from  those  noisy  signals. 

A  Kalman  filter  uses  all  indirect  measurements  of  data,  plus  prior  knowledge  of 
system  and  covariance  information  of  both  the  state  variable  and  indirect  measurements, 
to  estimate  the  desired  state  variable  in  a  manner  that  minimizes  the  error  statistically. 

From  a  Bayesian  viewpoint  [17,  20,  22],  the  Kalman  filter  propagates  the 
conditional  probability  density  of  the  desired  variable,  conditioned  on  actual  data  coming 
from  measuring  devices.  The  conditional  probability  density  propagation  depends  on 
assumptions  that  linear  models  can  represent  the  system  and  that  the  system  and 
measurement  noises  are  white  and  follow  a  Gaussian  (Normal)  distribution. 
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Figure  2.  A  Typical  Kalman  Filter  Application  [16:5] 


2.4.2  System  Models 

To  apply  the  Kalman  filter,  we  generally  need  [8]: 

1)  A  dynamic  model  that  indicates  how  the  state  vector  changes  over  time. 

2)  A  measurement  model  that  relates  observations  to  the  state  vector  and  updates 
estimates  based  on  these  observations. 

3)  An  estimate  of  the  variance  of  the  dynamic  noise  covariance. 

4)  An  estimate  of  measurement  (or  observation)  noise  covariance. 

In  some  applications  we  don’t  have  estimates  of  both  the  dynamic  and 
measurement  noise  covariances.  For  a  series  of  scalar  values,  we  cannot  simultaneously 
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estimate  both,  since  they  are  indeterminate.  Instead  we  may  directly  estimate  the  ratio, 
which  is  the  Kalman  filter  gain.  The  Kalman  filter  gain  is  presented  in  Section  2.4.3. 

Before  considering  Kalman  filter  equations,  we  present  underlying  stochastic 
equations:  linear  stochastic  difference  equations,  and  associated  dynamic  models  and 
measurement  models.  Since  we  use  discrete,  time-invariant  Kalman  filters,  we  limit  the 
discussion  to  discrete-time  system  models  in  this  section  as  well. 


2.4.2. 1  Linear  Stochastic  Difference  Equation 

The  general  form  of  the  linear  stochastic  difference  equation  is 
X(?n  )  =  <&(tntn_x  )X(tn_l  )  +  Bd  )u(tn_x )  +  Gd  (t„_,  )wd  (tn_x )  (1) 

where 

tn  is  the  n,h  time  interval, 

the  indices  d  shows  that  they  are  discrete, 

X  is  a  vector  of  system  states, 

0  is  transition  matrix  of  states  from  one  time  to  next  time, 

Bd  is  discrete  system  control  input  matrix, 
u  is  discrete  deterministic  control  input  matrix, 

Gd  is  discrete  dynamic  noise  input  matrix, 

wd  is  zero-mean  white  Gaussian  noise  vector  with  covariance  Qd  (tn_x )  [16]. 
The  mean  and  covariance  of  the  X(tn)  process  defined  by  Equation  1  are 


p„  (».)=«•(».,  V.  )P»  *»’■  ) + o , )Q,  (t,_,  )G/  «... ) 
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2A.2.2  Dynamics  Model 

The  discrete-time  version  of  the  dynamics  equation  can  be  obtained  by  integrating 
the  continuous-time  dynamic  equations  over  time.  So  the  discrete-time  dynamics 
equation  without  control  inputs, 

X0„ )  =  <&{tntn_x  )X(tn_l  )  +  Gdwd  (tn_ j ) 

where 

wA-i)  =  f  r)G(T)dfi(T)  ,  and 

Jtn~  1 

/?is  dynamic  driving  noise,  which  can  be  modeled  as  Brownian  motion  [16]. 

The  transition  matrix  relates  the  state  vector  X(r„_, )  at  one  time  to  a  state 

vector  at  the  next  time  index  X(tn ) .  w(tn_, )  is  normally  distributed  with  mean  of  zero  and 
have  the  properties: 

£K<X,-i)]  =  o 

(C.  W  (C, )]  =  £_  ®(tnT)G(T)Q(T)GT  ( r)<t>T(tnT)dr  =  Qd  (r„_, )  (2) 

E[wd  (t,  )w/(tj )]  =  0  for  t,  *  tj 
Discrete-time  dynamics  equation  can  be  also  represented  as 

X(X>)  =  ®(!ntn_ ,  )X(?n_, )  +  Gd  (r„_,  )wrf  (r„_, )  (3) 

If  we  assume  an  equally  spaced  and  time-invariant  system  model,  the  transition 
matrix  is  constant,  that  is  ®(tn_Ltn)  =  d> .  The  dynamic  noise  wd(tn_1)  is  assumed  to  have 

constant  variance,  Qd  (tn_x )  =  Qd ,  and  also  noise  input  matrix  G(tn  , )  is  also  assumed  to 
be  constant. 
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The  time-invariant  state  dynamics  model  can  be  represented  as 

X(?n)  =  $X(v1)  +  GW(f„_1)  (4) 

We  use  this  dynamic  model  to  represent  ARMA  models  later  in  Chapter  3. 

2.4.2.3  Measurement  Model 

The  modeled  measurement  or  observation  vectors,  Z(t„)  ,  are  linearly  related  to 
the  state,  but  observed  with  measurement  errors: 

Z(0  =  H(fJX(0  +  v(0  (5) 

where 

Z  is  the  measurement  or  observation, 

H  is  the  measurement  or  observation  matrix, 

X  is  the  vector  of  system  states, 
v  is  the  measurement  noise. 

The  measurement  noise  is  assumed  to  be  normally-distributed  white  noise  with 
zero  mean,  such  that 

E[v  (/„)]  =  0 

E[v(tn)vT(tn))  =  R(tn),  and 
E[v(tt )v  T (tj )]  =  0 for  t,*tj 

The  covariance  of  the  time-invariant  measurement  white  noise  is  R.  The  dynamics 
noise  and  measurement  noise  are  assumed  to  be  uncorrelated.  The  measurement  matrix 
H(tn )  and  measurement  noise  process  R(tn )  are  assumed  not  to  change  with  time. 
Therefore,  they  are  replaced  with  H  and  R  respectively.  The  time-invariant  measurement 
model  can  be  expressed  as 
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Z(0  =  HX(0  +  v(0 


(7) 


We  employ  this  model  assuming  the  ARMA  data  are  observed  with  additive 
observation  noise. 

2.4.3  Discrete  Space-Time-Invariant  Kalman  Filter  Algorithm 

In  this  section,  we  present  the  two  stages  of  Kalman  filtering:  the  propagation 
stage  and  the  measurement  update  stage.  First  let  us  define  some  terms.  Discrete  means 
that  two  stages  occur  only  at  set  intervals,  as  opposed  to  continuous.  Time-invariant 
implies  that  the  matrices  <1 >,G,  and  H  do  not  change  throughout  the  process  (stationary 
process).  The  covariance  matrices  of  the  noise  terms  Qd  and  R  are  constant,  assuming 
the  process  is  stationary. 

At  the  time  of  an  observation,  two  state  estimates  exist.  These  state  estimates  are 
indicated  by  a  hat  over  the  state  vector.  The  estimate  prior  to  a  measurement  update  is 

/V  _ 

labeled  with  a  superscript  minus  sign,  \(tn  ) .  Similarly,  the  state  estimate  after 

'v  , 

measurement  update  is  labeled  with  a  superscript  plus  sign,  X  (tn  )  .  The  associated 

state  covariance  matrices  are  P (t~  )  and  P  (t+  )  ,  respectively. 

These  state  estimates  are  the  conditional  distribution  of  the  state  vector.  With  the 
Gaussian  (normally  distributed)  assumption,  the  state  estimates  and  associated  covariance 
matrices  specify  the  conditional  distribution  of  the  true  state  (in  a  Bayesian  sense). 

2.4.3. 1  Propagation  Stage 

The  Kalman  filter  iterates  between  propagations  through  time  and  updates  for  each 
available  measurement.  The  Kalman  filter  gain  weights  the  dynamic  model  versus  the 
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measurement  information  when  updating  estimates.  The  state  estimate,  which  is  the  mean 
of  the  conditional  state  distribution,  propagates  through  time  with 

X(0  =  )X(tM_j+  )  +  Gdwd  (t„_j )  (8) 

where 

X  is  the  estimated  state  vector, 

<t>  is  the  transition  matrix, 

Gd  is  the  dynamic  noise  input  matrix. 

The  associated  growth  in  the  covariance  matrix  determined  by 

P(0  =  «*(f+-*)®r  +  GdQdGdT  (9) 

where 

P  is  covariance  matrix  of  state  estimates, 

Gd  is  the  dynamic  noise,  wd  (tn ) ,  input  matrix,  and 

Qd  is  the  covariance  matrix  of  discrete  dynamic  driving  noise. 

2A.3.2  Measurement  Update  Stage 

The  measurement  stage  combines  two  sets  of  information:  a  state  estimate  based 
on  the  dynamics  model  in  the  propagation  stage  and  a  correction  to  this  estimate  based  on 
an  actual  measurements  and  the  measurement  model.  For  known  dynamics  noise 
covariance  matrix  Q  and  known  measurement  noise  matrix  R,  the  state  covariance 
matrices  are  independent  of  the  actual  observations.  The  Kalman  filter  provides 
weighting  between  the  two  sets  of  information  about  the  state.  Kalman  filter  gain  is 
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K(0  =  P(f„-)HT[HP(C)HT  +  R]-1 


(10) 


where 

H  is  measurement  or  observation  matrix. 

After  measurement  at  time  tn ,  the  updated  state  estimate  is 

X(C  )  =  X(C)  +  K  ( tn  )[z (tn )  -  HX(<; )] 

where 

Zn  is  actual  measurement,  and 

HX  (tn  )  is  the  prediction  of  measurement  based  on  assumed  measurement 

model.  If  we  define  the  residual  vector  as  the  measurement  minus  the  expected 
measurement 

r(l.)  =  «(i,)-HX(i;)  (11) 

we  can  rewrite  the  update  equation  as 

£(C)  =  £((„-) +  K(0r(0  (12) 

The  update  of  the  covariance  relationship  is  calculated  with 

P(r;)  =  P(0-K(r„)HP(r„-)  (13) 

Again,  for  known  dynamics  noise  covariance  matrix  Q  and  known  measurement  noise 
matrix  R,  the  state  covariance  matrixes  are  independent  of  the  actual  observations. 
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Figure  3.  Kalman  Filter  Process  [8] 


In  many  engineering  applications,  these  covariances  matrices  are  pre-computed 
and  stored.  A  single  propagation  and  update  stages  of  a  Kalman  filter  process  is 
illustrated  in  Figure  3. 

2.5  Multiple  Model  Adaptive  Estimation  Algorithm 

The  theoretical  development  of  MMAE  is  given  in  Section  2.1.  In  this  section,  we 
present  the  algorithms  of  MMAE.  The  theory  and  notation  presented  in  this  section 
closely  follow  the  development  by  Maybeck  [16, 17]. 
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Let  a  denote  a  vector  of  parameters  to  be  estimated  (by  MMAE),  which  affect  any 
or  all  of  ^>,Bd,H,Qd,  and  R .  The  continuous  range  of  values  for  a  is  discretized  into  K 

representative  sets  of  values  (a,,a2, . ,ak ),  where  each  ak  corresponds  to  model  with  kth 

Kalman  filter.  The  Equations  1,  2,  5,  and  6  with  the  A:'*  filter  model  can  be  given  as 

X*  (O  =  *  *  ( V- A  (*-l  )  +  G<lk  ('«-!  (* -1 ) 

Z*(0  =  H *(0x*(0+v*(0  with 
)wrf/  )]  =  Qdk  (L_, ) 

^v4(f>4r  (*,,)]  =  MO 

Based  on  this  model,  the  Kalman  Filter  propagation  and  measurement  update 
Equations  8, 9, 10, 11, 12,  and  13  are  also  given  with  addition  of  subscript  k  on  all 
variables 

X*(0  =  *kit^)Xk(t^)  +  G*wdk{t,J 
p*  (O  =  (*Vi)<I>/  +GdkQdkGdkT 

rk(K)  =  Z  *  (^n  )  “  H  k  X  k  (t  n  ) 

Kt(g  =  pi(OH;[H4pt(t;)HiT+R4r1 

A4(r.)  =  H4p4(f;)H4T+R  t  ,  which  is  variance  matrix  of  measurements 

X4(0  =  X*(Q  +  K*(0r*(0 

Pk(t+n)  =  Pk(rn)-Kk(tn)HkPk(0 
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Each  Xk  (tn+ )  is  then  calculated  using  a  Kalman  filter  with  the  associated 
ak  parameters.  The  probabilities  are  calculated  using  their  residuals 
rk  (tn )  =  z*  (t„ ) _  X*  (t~ )  and  their  covariance  matrices  Ak(tn)=HkPk(t~)HkJ  +Rk  . 

Since  these  are  jointly  normally  distributed,  the  likelihood  of  each  filter  is  calculated  by 


f Z(t„)\a,Z(tnA)(Zn 


ak,Zn_,)  =  {2n) 2 1 A, 1 2  expjy  rkT  (tn)Ak\  (r„)j 


(14) 


where  S  is  the  dimension  of  the  measurement  vector  Z„,  S=1  in  this  application. 

Let  the  probability  that  a  assumes  the  discrete  value  ak ,  conditioned  on  the 

observed  measurement  history  to  time  tn  be 

Pk^n)  =  Probla  =  a*|z(0  =  ZJ 

Then  pk  (tk)  can  be  evaluated  recursively  for  all  k  via  the  iteration 


Pk  (0  = 


f  Z(fn)|fl,Z(fn_i)  (Zn  \ak  ^n-\)Pk  (tn- 1  ) 

f Z(r„)|a,Z(r„_,)  (Zn  ak  ’  Z„-1  )pk  (K- 1 ) 


(15) 


The  largest  likelihood,  calculated  in  Equation  14,  based  on  a  residual  in 
consonance  with  its  covariance  gives  the  best  set  of  parameters  with  a  high  probability 
from  Equation  15. 

The  probabilities  at  each  measurement  time,  tn,  for  n  =  1,...,  N,  must  sum  to  one: 

*=1 

This  normalization  limits  the  conditional  probabilities  to  only  the  K  discrete  parameter 
combinations  used  in  the  filters. 
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The  Bayesian  minimum  mean  squared  error  state  estimate,  the  sum  of  the  K 
discrete  state  estimates  weighted  by  their  associated  probabilities,  is  the  MMAE 
estimated  state  vector: 

X*„«(C)  =  EX,  «/)/>,('„)  (16) 

k= 1 

The  covariance  matrix  for  the  updated  MMAE  state  vector  XMMAE  (tn+)  is 

k= 1 


and  the  conditional  mean  for  the  unknown  vector  of  system  parameters  a  at  time  tn  is 
a(tn  )  =  E{a-  a(tn )  Z(tn  )  =  Zn}  =  ^akpk(tn) 

k=l 
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CHAPTER  III:  METHODOLOGY 


3.1  Overview 

Having  already  presented  the  detailed  general  MMAE  algorithm  in  Section  2.5, 
this  chapter  presents  its  application  to  time  series  analysis.  Since  the  MMAE  algorithm  is 
taken  from  engineering,  and  this  is  the  first  attempt  to  apply  it  to  a  time  series,  we  explain 
it  from  an  operational  research  point  of  view.  We  present  state  space  models  in  Section 

3.2  and  Kalman  filter  Section  3.3.  Finally  Section  3.4  presents  the  MMAE  forecasting 
technique 

3.2.  State  Space  Models 

So  far,  we  have  assumed  that  the  observed  time  series  measurements  are  exact. 
Therefore,  the  classical  and  Box- Jenkins  techniques  assume  that  at  time  t  the  time  series 
actually  had  the  exact  value  of  X,.  Box,  Jenkins  and  Reinsel  [3]  note  that  we  may  employ 
models  that  assume  superimposed  or  additive  white  noise.  Equation  7  can  express  the 
corresponding  system  observation,  where  X,  is  a  multivariate  time  series  with  the 
components  called  state  variables,  H  is  the  mapping  matrix,  and  v  is  white  noise. 

We  assume  that  state  X,  depends  on  the  previous  state  as  in  Equation  8,  from 
which  we  delete  the  time  indices  of  transition  matrix,  due  to  time  invariance.  This 
equation  is  called  the  state  equation.  Both  the  state  equation  and  the  measurement 
equation  form  a  state  space  model.  In  time  series  perspective,  X(tn )  represents  a  time 

series  that  follows  an  ARMA  structure.  In  the  next  section  we  show  how  to  represent 
ARMA  models  in  state  space  form. 
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3.3  Kalman  Filters  Method 


Even  if  we  use  an  ARMA  dynamics  model  to  represent  the  state  of  the  system,  X,, 
the  presence  of  observation  noise  induces  two  major  changes  to  the  forecasting  process. 
The  first  change  results  from  a  lack  of  precise  knowledge  regarding  the  state  of  the 
system  at  a  given  time.  We  discuss  the  probability  distribution  of  the  true  unobserved 
state.  The  second  change,  we  must  choose  how  to  treat  this  true,  but  unknown,  state.  We 
take  a  Bayesian  approach  (imbedded  in  Kalman  Filter)  that  this  unknown  state  is  a 
random  variable.  Therefore,  we  begin  with  a  prior  distribution  of  the  state,  and  update 
this  distribution  to  a  posterior  distribution  based  on  the  information  provided  by  each 
observation. 

The  Kalman  filter  recursive  equations  enable  easy  calculation  of  one-step 
forecasts,  provided  a  model  is  written  in  state  space  form.  The  details  of  Kalman  Filter 
equations  are  given  in  Chapter  II.  In  this  chapter,  we  give  the  corresponding  matrices  and 
apply  them  to  time  series  data. 

3.3.1  The  State  Dynamics  Model 

The  state  dynamics  model  in  this  research  is  ARMA  (2,2)  model.  We  give  the 
state  space  formulation  of  an  ARMA  (2,2)  model  below.  The  general  form  of  mean- 
corrected  ARMA  (2,2)  model  is  normally  written  as 

X,  =  *XM  +  #,_2  +£,~  #,£,-!  -  62£,-2 

where  X,  =  X,  -  X ,  with  X  =  —  V  Xf  .  Then  we  can  write  ARMA  (2,2)  model  in  state 

n  i=i 

space  form,  X(tn)  =  d)X(t,1_1 )  +  Gw(tn_{ ) ,  with 
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where  w(rn_j )  is  normally  distributed  with  zero  mean  and  variance  of  dynamic  noise 
Q  =  ct2£. 


3.3.2  The  Measurement  Model 

The  measurement  or  observation  vector,  Z  ( tn ),  Equation  7,  is  linearly  related  to 
the  state,  but  observed  with  measurement  noise,  where  the  time-invariant  measurement 
white  noise,  v(tn) ,  is  normally  distributed  with  zero  mean  and  covariance  R.  The 

dynamic  noise,  Q,  and  measurement  noise,  R,  are  assumed  to  be  uncorrelated.  In  our 
application  each  noises  and  their  variances  are  scalars. 

If  the  time  series  is  the  first  element  in  the  state  vector  of  dimension  m,  the 
measurement  vector  H  has  m  terms,  [1,  0,  0,  ...,0]  and  the  measurement  noise  is  a  single 
value.  In  our  application  m  is  equal  to  4. 

When  the  measurement  model  is  used  with  an  ARMA  dynamics  model,  the 
modeled  process  is  more  general  than  a  traditional  ARMA.  It  is  in  fact  an  ARMA  time 
series  observed  with  measurement  noise.  The  dynamics  noise  affects  future  observations, 
while  the  measurement  noise  only  affects  the  observation  at  that  time.  In  fitting  an 
ARMA  observed  with  measurement  noise,  the  data  series  are  used  as  the  realized 

observations,  Zn  ■ 
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3.3.3  Kalman  Filter  Forecasting 

In  applying  a  single  Kalman  filter,  we  go  through  the  same  four  steps  as  the  Box- 
Jenkins  methodology:  identification,  estimation,  diagnostics,  and  forecasting. 

3.3.3.1  Model  Identification 

For  ARMA  dynamics  models,  model  identification  is  the  same  as  the  Box -Jenkins 
methodology.  The  details  exceed  the  scope  of  this  thesis,  but  basically  model 
identification  is  done  by  looking  at  the  patterns  of  both  autocorrelation  and  partial 
autocorrelation  functions  of  data.  As  mentioned  before,  we  accept  ARMA  (2,2)  as  our 
dynamics  model  [2,15].  In  fact,  Box-Jenkins  [3]  and  Makridas  et.al.  [15]  do  not  suggest 
using  more  parameterized  model  than  an  ARIMA  (2,2,2),  which  covers  a  tremendous 
range  of  practical  forecasting  situations.  Moreover  the  parsimony  principle  [3]  tells  us  to 
use  the  minimum  possible  amount  of  parameters. 

The  selected  ARMA  (2,2)  model  is  used  for  the  dynamics  equation  in  the  Kalman 
filter.  Since  parameters  may  be  estimated  to  be  zero,  the  degenerative  ARMA  models  are 
also  possible. 

3.3.3.2  Model  Estimation 

As  in  the  ARIMA  estimation,  we  conduct  a  nonlinear  search  of  the  parameter 
space.  The  Kalman  filter  parameter  space  has  one  more  parameter  than  the  underlying 
Box-Jenkins  ARMA  model.  As  in  the  ARMA  model,  we  have  <j>i,  <j>2, ...,  <J)P,  0i,  02, ..., 

Qq,  as  appropriate,  and  Q  =  (dynamic  noise  variance  in  Box  Jenkins  and  Kalman 
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filter  notation).  The  additional  parameter  is  either  the  first  term  of  the  Kalman  filter  gain, 
k,  or  the  measurement  noise  variance  R. 

We  may  determine  the  Kalman  filter  gain  using  the  set  of  ARMA  coefficients  to 
determine  the  transition  matrix  d>  along  with  the  dynamic  noise  variance  Q  and  the 

measurement  noise  variance  R.  Begin  with  any  initial  estimate  of  the  P(/0 ) ,  we  may 

iterate  through  the  propagation  and  update  equations.  Equations  9, 10,  and  13,  until  the 
Kalman  filter  gain  becomes  constant. 

Notice  that  the  methodology  depends  upon  both  the  dynamics  noise  shock 
variance  Q  and  the  measurement  noise  variance  R.  While  each  noise  is  assumed  to  be 
independent,  their  impact  on  the  Kalman  filter  gain  is  highly  correlated.  Maybeck  [21] 
states  that  ratio  of  the  magnitude  of  the  eigenvalues  determines  the  magnitude  of  the 
Kalman  filter  gain.  In  our  application,  since  the  noise  and  their  variances  are  scalars,  the 
ratio  of  Q  and  R  (for  any  given  set  of  model  coefficients)  determines  the  steady  state 


Kalman  filter  gain,  which  is  simply  K  = 


Q 


Q  +  R 


Therefore,  we  fix  Q=1  and  treat  R  as  a 


multiple  of  Q.  The  true  Q  may  be  estimated  from  the  residuals,  if  desired. 

We  apply  the  ARMA  dynamics  model  by  iteratively  propagating  and  updating  for 
all  the  available  observations.  We  propagate  the  state  estimates  using  Equation  8.  We 
update  the  state  based  on  the  next  observation  and  a  steady-state  gain  with  Equations  1 1 
and  12.  The  residuals  (equivalently  errors)  are  used  to  evaluate  each  set  of  coefficients. 

We  conduct  a  nonlinear  search  to  determine  the  model  coefficients  and  R.  The  set 
of  model  coefficients,  which  maximize  the  maximum  likelihood  function  and 
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equivalently  the  conditional  joint  distribution  of  observations,  shown  in  Equation  17 
[3:276-277],  is  selected  as  the  best  fitting  model. 

p(z  |  (p,0,Q)  =  ft'  (2^ a,))"' 7 2 exp 

._  -i  /  2  r 

=  exp  -  xr«,)j 

i=l  LZ/11V,/  1=1 

We  calculate  the  residuals  with  Kalman  filter  equations.  Since  only  the  ratio  of  Q  and  R 
is  important,  we  take  Q  =  1  (arbitrarily,  we  can  take  any  value)  and  R  any  starting  value 
greater  than  zero. 

During  the  nonlinear  search  phase,  to  eliminate  the  initialization  bias,  we  first  find 
the  initial  state  estimates  by  backward  forecasting  or  backforecasting.  To  begin  the 
Kalman  filter  recursion  for  backforecasting  and  forecasting  we  need  an  initial  state 
estimate  and  an  initial  covariance  matrix.  We  use  the  expected  mean  corrected  state, 
which  is  the  zero  vector,  as  the  initial  state  estimate.  We  may  select  P(t0) to  be  very 

large,  to  represent  a  non-informative  prior.  If  we  desire  the  exact  maximum  likelihood 
estimates  (MLE),  for  an  ARMA  (2,2)  we  fix  R= 0,  which  treats  the  measurements  as 
exact,  and  solve  for  the  expected  covariance  based  on  our  state  definition.  We  may  use 
the  general  linear  filter  representation  to  determine  this  covariance  matrix.  Recall  the 
general  linear  filter  is 

f  I  =  X  t  ~  M  =  ^  W  k  £t-k 

k  =  0 

Since  the  errors  are  uncorrelated 

Var(*()=  £  w\  Var(£-t.J  =  cr2  £  \f/2k  ,  and 

k=0  k-0 
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Cov(x<,  Xt+k)=(T2  E  ViVi+k 

i~0 

For  an  AR(2)  process,  the  weights  may  be  calculated  with  the  recursive  equation 
y/  =  (f)  1  \f/  .  j  +  (/>2  ty  j  2  start*n§  with  ¥  o  =  1  and  if/  _x  =  0  .  The  coefficient  of  the 

moving  average  terms  may  be  added  directly  to  the  appropriate  filter  coefficient  later. 

After  calculating  all  weights  for  an  ARMA(2,2)  model,  with  the  above  recursion, 

the  first  two  weights  are  recalculated  with  If/  x  =  0 1  +  0X  and  lf/2  —  01\f/1  +  02  +  62. 

Thus,  we  may  determine  the  initial  state  covariance  matrix  as,  for  an  ARMA  (2,2) 
dynamics  model  with  state  definition  of 

5' 

x(rj= 

e, 

J,- 1 

with 

Var{X,}  cov  { X, ,  X,_, }  cov{X„e,}  covlZ,,^.,} 

cov{Z„XM}  Var{X,_x}  cov  {*,_„£,}  cov{X,_pfM} 
cov{A,,f,}  cov{X, Var{st}  cov{f, 

covJX,,^}  CO  v{X,.,,  £•,.,}  cov  {£•,,£•,.,}  Var{e,_x} 

Var{X, }  cov{X„XM}  Q  ysx2Q 

cov{X(,Z,_j}  Var{X,_x )  0  y2Q 

Q  0  Q  0 

Vx2Q  VxQ  0  Q 
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£  V‘2  £  1  Vx 


i=0 


i=0 


£  £  Vi 


i  =  0 


0 


0  ^ 

1  0 
0  1 


(18) 


Therefore,  with  the  Kalman  filter  equations,  we  may  obtain  the  MLE  values  for 
ARM  A  models  [10]. 


3.3.3.3  Forecasting 

After  the  model  parameters  are  selected,  beginning  with  initial  state  estimates 
equal  to  zero  vector,  and  an  initial  covariance  matrix  P(t0 )  as  calculated  in  Equation  18, 
we  process  the  available  data  through  the  Kalman  filter  using  recursive  equations.  When 
we  complete  the  last  datum,  the  last  state  estimate  is  conditioned  on  all  prior 
observations.  Therefore  future  forecasts  are  all  conditional  estimates. 

The  one-step  ahead  forecast  is  available  directly  from  the  propagation  equation: 

x(c,)  =  ox(0 

Predictions,  /-steps  into  the  future,  require  propagation  of  these  estimates,  without 
the  update  stage  with 

x(<;«  I  z(<„ )  =  z„ ) = <j>'  x(C) 

All  Kalman  filter  state  estimates  are  conditional.  In  the  prior  equation,  the 
conditioning  is  explicitly  included  since  it  does  not  follow  the  convention  of  conditioning 
all  observations  prior  to  the  time  of  the  state  estimate. 

While  missing  data  is  difficult  to  analyze  with  classical  and  traditional  Box- 
Jenkins  techniques,  the  Kalman  filter  can  easily  handle  missing  observations.  The  best 
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state  estimate  is  used  in  place  of  the  missing  observations.  Therefore,  we  propagate 
without  an  update  for  missing  observations. 

3.3.3.4  Prediction  Intervals 

Based  on  the  estimated  dynamics  noise  variance  and  the  measurement  noise 
variance  (or  equivalently  the  Kalman  filter  gain),  we  may  calculate  all  Kalman  filter  state 
covariance  matrices.  If  we  assume  the  states  follow  a  Normal  distribution,  we  may  make 
a  Bayesian  probability  interval  for  future  forecasts.  The  mean  is  the  state  estimate  based 
on  the  data  available  and  the  variance  is  the  corresponding  covariance  matrix. 

We  use  empirical  prediction  errors  for  the  available  data  set  to  estimate  prediction 
variance.  We  process  the  data  through  Kalman  filter  once  more,  this  time  making  1-step, 
3-step  and  5-step  ahead  forecasts  within  the  data.  All  prediction  errors  are  then  collected 
and  the  MSE  (mean  square  errors)  of  forecasts  found  for  every  step.  Using  these  MSE’s 
and  standard  prediction  interval  equation  we  find  the  prediction  intervals. 

3.4  Multiple  Model  Adaptive  Estimation  Forecast 

The  detailed  algorithm  is  given  in  Chapter  2.  In  this  chapter  we  present  the  time 
series  forecasting  application  of  MMAE. 

MMAE  applies  a  bank  of  Kalman  filters  to  process  the  data  simultaneously.  Its 
major  advantages  are:  first  it  combines  the  identification  and  estimation  step,  second  it 
probabilistically  weights  the  models,  and  gives  more  probability  to  the  more  likely 
model. 

We  begin  MMAE  by  selecting  discrete  set  of  coefficients  that  cover  the  feasible 
parameter  region,  rather  than  searching  for  a  single  set  of  parameters  that  best  fits  the 
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data.  For  the  AR  (2)  and  MA  (2)  portions  of  ARMA  (2,2)  dynamics  model,  the 
coefficients  are  selected  every  0.1  interval  for  <J>i,  (J>2  0i,  and  02  in  their  invertible 
triangle,  Figure  4.  The  bank  of  Kalman  filters  consists  of  a  single  Kalman  filter  at  each 
combination  of  these  parameters. 
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Figure  4.  Parameter  Selection  from  Invertible  Triangle  for  AR  and  MA  Models 


Each  Kalman  filter  processes  data  independently  of  each  other,  as  discussed 
previously.  Each  filter  starts  conditionally  with  an  initial  state  estimate  of  zero  vector, 
and  initial  covariance  vector  P  (to).  All  required  equations  are  given  in  Chapter  2. 

The  MMAE  algorithm  calculates  the  probabilities  associated  with  each  of  the  K 
Kalman  filters.  Based  on  the  assumption  of  zero  mean  and  the  estimated  residual 
variance,  the  normal  probability  density  function  for  the  nth  measurement,  zn ,  conditioned 
on  the  kh  filter's  vector  of  parameters,  a*  and  the  prior  measurement  history,  Z„_i  ,  is 
calculated  with  Equation  14. 
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The  probability  for  the  kth  filter  having  the  "correct"  parameters  conditioned  on 
the  measurement  history  through  time  t„  is  calculated  with  Equation  15. 

The  initial  or  a  priori  probabilities  account  for  information  available  about  the 
likelihood  of  particular  filter  combinations  before  the  measurement  data  are  processed.  If 
no  information  is  available,  the  a  priori  probabilities  should  all  be  equal;  pk  (to)  =  1  IK  for 
k  =  1,...,  K.  In  addition,  if  any  of  the  filter  probabilities  become  zero,  that  filter's 
probabilities  remain  zero.  To  prevent  prematurely  discarding  potentially  viable  filter 
parameters,  practitioners  commonly  apply  a  heuristic;  if  any  of  the  filter  probabilities 
decreases  below  a  very  small  lower  bound,  such  as  0.0001,  the  heuristic  artificially 
increases  that  filter’s  probability  to  the  lower  bound.  The  filter  probabilities  that  result 
after  the  last  datum  are  not  adjusted  with  this  heuristic.  The  final  filter  probabilities 
represent  the  likelihood  of  each  combination  of  model  parameters  conditioned  on  the 
available  measurement  history,  Z^. 

After  determining  the  probabilities  based  on  all  available  observations,  we  have  a 
conditional  (Bayesian)  probability  of  each  filter’s  parameters  being  correct.  We  use  these 
probabilities  to  calculate  a  conditional  estimate  by  probabilistically  weighting  each 
filter’s  forecast.  The  mean  estimate  of  next-step  forecast  conditioned  on  the  available 
measurement  history,  Z n,  is  calculated  with  Equation  16. 

All  Z-step  ahead  forecasts  for  each  filter  can  be  made  the  same  way  the  single 
Kalman  filter  does,  where  for  each  filter, 

*((';jz(0  =  z,)  =  ®'x(0 

Using  the  latest  calculated  filter  probabilities  we  can  weight  each  Z-step  forecast 
and  find  MMAE  forecasts  using 
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K 

^  MMAE  (tn+l  )  =  (*«+/  )P*  (*n  ) 

*=1 

Using  covariance  propagation  equation  without  the  update  stage,  we  calculate  all 
variances  for  each  /-step  forecast  of  each  Kalman  filter  as 

P*(^')  =  ®P*U+«Hi-i)®r+Gj3Gr 

Since  we  stop  calculating  covariance  matrices  after  the  Kalman  filter  gain 
becomes  steady,  we  use  the  last  calculated  P(t+«  )  of  each  filter  for  the  first  step 
covariances  in  the  above  equations. 

Based  on  calculated  covariances  we  find  the  prediction  interval  of  every  single 
forecast  for  every  filter,  then  weight  them  with  the  filter  probabilities  for  every  step.  The 
weighted  sum  of  the  lower  bound  and  weighted  sum  of  the  upper  bounds  of  prediction 
intervals  give  the  final  prediction  intervals  for  MMAE  forecasts. 
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CHAPTER  IV:  ANALYSIS  AND  RESULTS 


4.1  Overview 

In  previous  chapters  we  proposed  a  new  method,  MMAE,  for  time  series 
forecasting.  Chapter  III  gives  the  detailed  algorithm  of  MMAE.  We  claim  that  MMAE  is 
at  least  as  good  as  other  well-known  classical  and  Box-Jenkins  models  in  literature.  In 
this  chapter  we  conduct  a  Monte  Carlo  simulation  to  support  our  claim,  and  compare  the 
results  of  MMAE  with  the  other  known  forecast  methods’  results.  Section  4.2  presents 
the  implementation  of  the  Monte  Carlo  simulation.  Section  4.3  shows  how  we  evaluate 
our  method,  Section  4.4  presents  which  statistical  measures  we  use  to  check  the  model 
accuracies  and  interprets  simulation  outputs,  and  finally  Chapter  4.5  gives  the  results  and 
findings. 

4.2  Design  of  Monte  Carlo  Simulation 

The  design  of  Monte  Carlo  study  should  be  carefully  optimized  to  cover  all 
parameters,  variables  or  criteria  affecting  the  model.  In  our  research  there  are  four  main 
areas  of  concern:  sample  size,  noise  variances,  data  generation  parameters  and 
techniques,  and  simulation  trials,  which  are  to  be  chosen  carefully  for  inferential 
purposes. 

As  Box-Jenkins  [3]  suggest  the  sample  size  at  least  100  to  get  satisfactory 
forecasts,  we  select  sample  sizes  to  be  at  least  equal  to  or  more  than  100.  Moreover  we 
try  8  different  sample  sizes  with  the  same  parameters  and  same  data,  to  measure  the 
effect  of  sample  sizes  in  MMAE  forecasting.  The  simulation  results  suggest  no  affect  of 
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sample  size  in  MMAE’s  forecasting.  We  use  different  sample  sizes,  usually  between  200 
and  400.  We  become  consistent,  and  use  the  same  sample  sizes  for  the  same  data 
generation  parameters  with  different  variances. 

To  measure  the  performance  of  MMAE  against  data  fluctuations,  we  generate 
data  with  two  different  variances,  big  and  small.  The  data  are  normally  generated  with  the 
ARMA(2,2)  equation  with  additive  noise  variance,  v(tn ) , 

X,  =  [l  +  ^XM  +  02X,_ 2  +  —  Q\^t-\  ~  +  vt 

where  dynamic  noise,  w{tn ) ,  from  Equation  4,  and  measurement  noise,  v(tn ) ,  from 
Equation  7  are  distributed  normally  with  mean  and  variances  as  shown  in  Table  1 . 


Table  1.  Dynamic  and  Measurement  Noise  Variances 


Dynamic  Noise 

Measurement  Noise 

Htn)=£t~N(  0,1) 

v(tn)~N( 0,0.1) 

MtB)=£,~N(  0,25), 

v(tn)~N(  0,1) 

Data  generation  parameters  are  chosen  to  provide  representation  throughout  the 
stationarity  and  invertibility  regions.  We  test  our  method  by  generating  many  different 
data  sets  using  different  combinations  of  AR  and  MA  parameters  from  their  parameter 
regions,  the  invertible  triangle,  as  shown  in  Figure  4  and  Figure  5.  We  also  generate  data 
with  pure  autoregressive  models,  AR(1)  and  AR(2),  pure  moving  average  models,  MA(1) 
and  MA(2),  and  more  parameterized  mixed  model,  ARMA(3,3),  to  test  the  robustness  of 
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MMAE.  The  number  of  cases,  time  series,  from  each  data  generation  technique  is  given 
in  Table  2.  Data  generation  parameters  and  related  simulation  results  are  given  in 
Appendix  A  and  B. 

Table  2.  Number  of  Cases  for  Each  Data  Generation  Model 


ARMA(2,2) 

AR(1) 

AR(2) 

MA(1) 

MA(2) 

ARMA(3,3) 

Total 

36 

10 

10 

7 

7 

6 

76 

We  split  the  data  set  into  two  parts,  an  initialization  set  and  a  test  set  or  hold  out 
set.  Most  of  the  data  are  used  as  initialization  set  to  estimate  any  parameters  and/or 
initialize  the  method.  Forecasts  are  made  for  the  test  set.  Since  the  test  set  is  not  used  in 
model  fitting,  these  forecasts  are  genuine  forecasts  made  without  using  the  values  of  the 
observations  for  these  times.  All  statistics  are  computed  for  the  errors  in  the  test  set  only. 

We  run  the  simulation  1000  times  for  each  parameter  set,  in  order  to  get  accurate 
results.  The  simulation  results  tables  are  the  averages  of  a  1000  replications.  Since  we 
have  76  cases,  our  evaluation  consists  of  76,000  data  series. 

4.3  Evaluation  Criteria 

The  evaluation  our  method  is  another  important  aspect  of  the  research.  We  can 
simply  fit  MMAE  to  data,  make  forecasts,  collect  the  statistics,  and  make  inferences  from 
those  statistics.  However  none  of  these  statistical  measures  give  a  good  basis  of 
comparison.  Does  an  MSE  of  6  or  SSE  of  56  indicate  a  good  or  bad  forecasting 
performance ? Another  basis  is  to  compare  MMAE  results  with  the  other  well-known 
sophisticated  methods’  results.  The  relative  comparison  of  results  obtained  from  MMAE 
and  other  methods  is  more  meaningful  and  useful  than  simply  computing  statistics  of  the 
MMAE.  So  we  decide  to  fit  other  methods  as  well  as  MMAE  to  data.  The  other  forecast 
methods  compared  with  MMAE  are: 

•  Naive 

•  Moving  Average  with  window  size  of  5 

•  Moving  Average  with  window  size  of  10 

•  Simple  Average 

•  Exponential  Smoothing 

•  Regression 

•  AR  (1)  ULS,  unconditional  least  square  model 

•  AR  (2)  ULS  model 
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•  MA  (1)  ULS  model 

•  MA  (2)  ULS  model 

•  ARM  A  (1,1)  ULS  model 

•  ARMA  (2,2)  model  with  ‘TRUE’  data  generating  parameters, 

•  ARMA  (2,2)  model  with  ‘Initial  Parameters  Estimates’  according  to  Box -Jenkins 

[3:202] 

•  ARMA  (2,2)  MLE,  maximum  likelihood  mode 

•  Single  Kalman  Filter  with  dynamic  model  of  ARMA  (2,2) 

Naive  method  simply  takes  the  most  recent  observation  as  forecast.  Moving 
average  methods  take  the  average  of  5  and  10  most  recent  observations  as  forecast. 

Simple  average  is  the  average  of  all  previous  observations.  Exponential  smoothing  is  the 
weighted  average  of  the  estimate  made  at  previous  time  and  the  observation  at  this  time. 
Regression  uses  the  standard  regression  methods,  by  regressing  the  observations  to  time. 
Other  methods  except  Kalman  filter  method  are  standard  Box -Jenkins  models.  Further 
details  about  the  methods  can  be  found  in  [3, 7,  15].  ULS  shows  that  all  necessary 
parameters  are  estimated  by  unconditional  least  square  method.  MLE  shows  the 
parameters  are  estimated  by  maximum  likelihood  estimates.  The  parameters  for  the  single 
Kalman  filter  method  are  also  calculated  using  the  maximum  likelihood  estimation 
procedure  [10]. 

In  this  research  we  compared  MMAE  estimates  with  many  other  methods’ 
estimates,  but  formally  we  take  ARMA  (2,2)  MLE  as  standard  as  proposed  by  time  series 
experts  [3,8,  10, 15].  We  give  the  analysis  of  results  in  the  Section  4.5  based  on  this 
comparison. 

4.4  Measuring  Forecast  Accuracy 

In  this  section  we  give  the  brief  description  of  statistics  calculated  to  compare  all 
methods. 
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All  of  the  statistics  are  based  on  fc-step-ahead  forecast  error  (made  at  time  t), 
which  is 

et+k  (0  =  Xl+k  —  Ft+k  ( t ) 

where  k  represents  all  future  times,  fc=l  to  n, 

X?+jt  is  &-step  ahead  observation. 


Ft+k  is  k- step  ahead  forecast. 

The  statistical  measures  we  use  in  comparison  are:  mean  error  (ME),  mean 
absolute  error  (MAE),  sum  of  square  errors  (SSE),  1-step  prediction  interval  and  its 
coverage  (PIW1MSE  and  PIC1MSE),  3-step  prediction  interval  and  its  coverage 
(PIW3MSE  and  PIC3MSE),  and  5-step  prediction  interval  and  its  coverage  (PIW5MSE 
and  PIC5MSE).  The  equations  for  these  statistics  for  a  single  series  with  n  forecast  are 
given  as 


1  V' 

ME=  2_J  en 
n  *=i 
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where  &k  =  Var[et+k  (?)]  for  all  t.  The  prediction  interval  computation  for  MMAE  is 

given  in  Chapter  III.  Prediction  interval  coverage  is  simply  the  number  of  times  the 
prediction  intervals  contain  real  data  across  5  forecasts  and  all  replications. 

Table  3  shows  an  example  of  a  Monte  Carlo  simulation  results;  where  true 
parameters  are  the  parameters  used  to  generate  data;  sample  size  is  the  number  of 
generated  data;  number  of  predictions  is  the  number  of  test  data  set  and  forecasts;  noise 
s td.  deviation  is  the  standard  deviation  of  randomly  generated  dynamic  noise;  error  std. 
deviation  is  the  standard  deviation  of  randomly  generated  measurement  noise,  which  is 
added  to  data.  The  percentages  at  the  bottom  of  table  give  the  percentages  of  replications 
that  MMAE  estimates  are  better  than  ARMA  (2,2)  MLE  estimates. 

All  the  results  in  Table  3  and  other  tables  in  appendices  are  averages  of  1000 
replication  for  that  data  generation  parameters  set.  In  this  case,  MMAE’s  forecasts  are 
better  and  more  accurate  than  MLE’s  forecasts  53.5%  and  55.5%  of  time  out  of  the  1000 
replications,  with  respect  to  statistics  MAE  and  SSE.  MMAE  gives  the  tightest  prediction 
intervals  with  high  coverages.  The  narrow  prediction  interval  express  the  forecast 
precision,  the  narrower  the  prediction  interval,  the  more  accurate  the  forecasts  are. 
Especially  the  MMAE’s  prediction  intervals  are  narrower  than  MLE’s  prediction 
intervals.  Further  conclusions  are  given  in  the  next  section. 
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Table  3.  Monte  Carlo  Simulation  Results 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  53.5  55.5 


4.5  Results  and  Findings 


All  Monte  Carlo  Simulation  Results  are  presented  in  Appendix  B. 

In  an  application  of  MMAE  for  a  time  series  we  recommend  using  all  possible 
Kalman  filters  within  the  parameter  space  of  AR  and  MA  parameters  in  their  invertible 
triangles  as  shown  in  Figure  4  and  Figure  5.  With  an  interval  of  0.1  and  ARMA(2,2) 
dynamic  model,  the  number  of  whole  possible  Kalman  filters  within  the  parameters  space 
is  about  104,976.  This  is  not  computationally  efficient  for  a  Monte  Carlo  Simulation, 
since  we  run  each  case  1000  times,  we  have  to  limit  ourselves  to  a  smaller  parameters 
space.  We  focus  on  parameters  around  the  optimal  parameters  of  single  Kalman  Filter 
ARMA(2,2)  forecasting.  Since  we  run  the  model  for  a  real  time  series  application  only 
once,  we  can  use  all  104,967  filters. 

We  examine  each  result  based  on  its  data  generation  model,  and  later  make  an 
overall  evaluation.  We  first  compare  the  averages  of  statistics  of  MLE  and  MMAE 
forecasting  techniques,  and  the  average  of  percentages  of  replications  MMAE  estimates 
better  than  MLE  estimates,  and  then  we  also  give  the  comparison  of  MMAE  with  the 
other  techniques  later. 

We  would  like  to  test  MMAE’s  capability  to  forecast  the  time  series  data  with  big 
and  small  variance.  We  run  the  same  parameterized  model  twice,  one  is  with  big  variance 
and  the  other  is  with  small  variance,  as  described  in  Section  4.3.  Here  we  give  the 
summary  of  those  results. 

Table  4  shows  that  MMAE  gives  4%  better  estimates  for  small  variance  than  big 
variance  for  ARMA(2,2)  generated  data,  both  compared  with  MLE.  These  percentages 
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alone  are  not  very  meaningful.  We  need  to  interpret  them  along  with  the  other  statistics 
averages  as  well. 


Table  4.  Percentages  of  MMAE  Estimates  Better  Than  MLE  Estimates  with 

ARMA(2,2)  Data 


MAE 

SSE 

Big  Variance 

40.5% 

39.9% 

Small  Variance 

44.4% 

43.5% 

Average 

42.5% 

41.7% 

From  Table  5,  ME’s  are  the  same,  MAE’s  and  SSE’s  are  very  close  to  each  other. 
We  also  see  that,  the  first  step  forecast  MMAE  gives  narrower  prediction  interval  with 
higher  coverage.  The  other  prediction  intervals  of  MMAE  are  also  narrower  with  slightly 
less  coverages  than  MLE. 


Table  5.  Averages  of  Statistics  for  ARMA(2,2)  Data  with  Big  Variance 


ME 

MAE 

SSE 

PIW1 

EMB 

PIW5 

PIC5 

MLE 

-0.01 

1.25 

23.48 

2.23 

0.93 

2.70 

0.91 

1.20 

13.86 

2.67 

0.89 

Although  the  percentages  for  the  small  variance  are  larger  than  big  variance,  the 
statistics  of  MMAE  for  small  variance  ARMA(2,2)  data  are  slightly  worse  than  MLE’s  as 
shown  in  Table  6.  If  we  evaluate  Table  5  and  Table  6  together,  compared  with  MLE’s 
statistics,  we  can  conclude  that  MMAE  forecasts  for  big  variance  data  are  better  than 
small  variance  data.  So  we  have  to  be  careful  while  interpreting  the  results.  The 
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percentages  do  not  necessarily  mean  better  forecast  and  do  not  reflect  the  magnitude  of 
the  errors  exactly.  The  counter  counts  the  percentages  without  looking  at  the  quantity  of 
the  error;  it  even  counts  for  very  small  differences  at  the  bigger  decimal  places.  It  is  seen 
that  some  of  the  statistics  of  small  variate  data  are  relatively  larger  than  those  of  big 
variate  data.  From  this  discussion  we  can  conclude  that  MMAE  forecasts  are  better  for 
more  variable  data  than  less  variable  data. 


Table  6.  Averages  of  Statistics  for  ARMA(2,2)  Data  with  Small  Variance 


ME 

MAE 

SSE 

PIW1 

PIC1 

PIW3 

PIC3 

PIW5 

PIC5 

MLE 

-0.11 

5.67 

262.43 

9.88 

0.88 

12.62 

0.92 

12.25 

0.90 

MMAE 

-0.14 

5.81 

291.91 

9.83 

0.87 

10.57 

0.86 

10.89 

0.84 

Table  7  shows  the  averages  of  Table  5  and  Table  6.  For  the  overall  average  of  all 
ARMA(2,2)  data  statistics,  MMAE’s  are  slightly  worse  than  than  MLE’s.  However 
prediction  interval  width  for  the  first  step  forecast  is  narrower  with  the  same  coverage. 
The  other  prediction  intervals  are  also  narrower,  but  have  4%  less  coverage.  As  a  result, 
with  averages  of  42.5%  and  41.7%  from  Table  4,  we  can  conclude  that  MMAE  forecast 
for  an  ARMA(2,2)  data  are  about  equal  to  MLE  forecasts. 


Table  7.  Averages  of  Statistics  for  all  ARMA(2,2)  Data 


ME 

MAE 

SSE 

PIW1 

PIC1 

PIW3 

PIC3 

PIW5 

PIC5 

MLE 

BiwiaJ 

3.46 

142.96 

6.06 

0.89 

7.82 

0.92 

7.47 

0.90 

MMAE 

BUM 

3.51 

152.88 

5.84 

0.89 

6.53 

0.87 

6.78 

0.86 
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In  Table  8,  MAE’s  and  ME’s  are  about  equal,  SSE  of  MMAE  is  worse.  The  first 
step  prediction  interval  of  MMAE  is  slightly  larger  with  higher  coverage,  and  the  other 
prediction  intervals  of  MMAE  are  narrower,  but  have  lower  coverages.  The  percentages 
are  very  close  to  50%.  We  can  say  that  MMAE  forecasts  are  not  better,  but  equal  to 
MLE’s  forecasts  for  AR(1)  data. 


Table  8.  Averages  of  Statistics  and  Percentages  for  AR(1)  Data 


ME 

PIC3 

PIW5 

PIC5 

MLE 

-0.10 

2.94 

99.70 

5.05 

MMAE 

-0.11 

3.12 

118.48 

5.12 

0.90 

6.08 

0.88 

6.33 

0.87 

Percentage 

48.3% 

47.4 

Table  9  shows  that  MMAE  forecast  for  AR(2)  data  are  not  good  compared  to 
MLE  for  ARMA(2,2).  All  statistics  except  first  step  prediction  interval  width  and 
coverage  are  worse  than  MLE’s  statistics. 


Table  9.  Averages  of  Statistics  and  Percentages  for  AR(2)  Data 


ME 

MAE 

SSE 

PIW1 

PIC1 

PIW3 

PIC3 

PIW5 

PIC5 

MLE 

EH! 

3.36 

119.20 

5.75 

0.86 

0.92 

7.68 

mm 

MMAE 

EQB 

3.98 

273.00 

6.48 

0.87 

| 

0.87 

7.09 

Percentage 

41.1% 

41.0% 

Although  the  percentages  are  less  than  50%,  Table  10,  MMAE  forecasts  are  about 
equal  to  MLE  forecasts.  Only  MMAE’s  MAE  and  SSE  are  slightly  worse  than  MLE’s. 
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All  the  other  better  statistics  and  narrower  prediction  intervals  and  higher  coverages  show 
that  MMAE  is  very  good  in  this  case. 


Table  10.  Averages  of  Statistics  and  Percentages  for  MA(1)  Data 


ME 

PIW1 

PIC1 

PIW5 

PIC5 

MLE 

82.29 

4.65 

5.46 

0.91 

MMAE 

2.63 

84.38 

5.00 

5.50 

0.92 

Percentage 

48.4% 

47.9% 

In  Table  11,  MMAE  forecasts  are  not  better  but  not  worse  either.  The  first  three 
statistics  seem  to  be  larger,  but  the  prediction  intervals  are  better  with  almost  the  same 
coverages. 


Table  11.  Averages  of  Statistics  and  Percentages  for  MA(2)  Data 


ME 

MAE 

PIW1 

PIC1 

PIC5 

-0.05 

3.37 

5.70 

0.91 

7.28 

0.92 

MMAE 

-0.11 

3.61 

0.91 

6.87 

0.87 

7.02 

Percentage 

44.6% 

45.5% 

For  ARMA  (3,3)  data,  Table  12,  MMAE’s  performance  is  not  as  good  as  MLE’s. 
MMAE’s  ME,  MAE  and  SSE  appear  very  close  to  MLE’s,  the  prediction  intervals  are 
narrower,  but  with  worse  coverages. 
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Table  12.  Averages  of  Statistics  and  Percentages  for  ARMA(3,3)  Data 


MAE 

PIW1 

PIC1 

PIW3 

PIC3 

PIW5 

PIC5 

MLE 

4.30 

5.89 

0.85 

9.21 

0.88 

9.66 

0.89 

MMAE 

0.03 

4.74 

277.67 

6.68 

0.84 

7.64 

0.79 

7.95 

0.86 

Percentage 

42.1% 

41.9% 

Tables  4-12  show  the  averages  of  simulation  results  based  on  each  data  generation 
model.  Since  the  pattern  of  data  of  any  time  series  can  be  any  of  the  above  models  and 
MMAE  does  not  require  the  data  to  have  any  specific  pattern,  we  think  that  if  we  take  the 
averages  of  the  above  models,  we  can  reach  a  more  general  conclusion  about  the 
goodness  of  MMAE’ s  forecasting  capabilities. 

Table  13  shows  the  averages  of  all  76  different  data  generation  models;  all  of 
them  have  different  variance  and  sample  sizes.  The  ME  and  MAE  of  both  methods  are 
almost  the  same,  SSE  of  MMAE  is  slightly  worse  than  MLE’s.  The  first  step  prediction 
interval  of  MMAE  is  only  2%wider  than  MLE’s  but  with  the  same  coverage  probability. 
The  PIW3MSE  of  MMAE  is  14%  narrower  than  MLE’s  with  only  4%  worse  coverage. 
The  PIW5MSE  of  MMAE  is  8%  narrower  than  MLE’s  with  only  3%  worse  coverage 
probability.  The  narrower  prediction  intervals  automatically  affects  the  prediction  interval 
coverages,  make  them  smaller.  Overall  averages  show  that  MMAE  forecasts  are  not 
better  than  MLE  forecasts,  but  they  are  not  worse  either,  they  are  almost  equal  to  each 
other. 


Table  13.  Overall  Averages  of  Statistics  and  Percentages  for  MMAE  and  MLE 


ME 

MAE 

SSE 

PIW1 

PIC1 

PIW3 

PIC3 

PIW5 

PIC5 

MLE 

3.07 

118.89 

5.21 

0.88 

6.87 

6.72 

0.90 

MMAE 

3.24 

5.33 

0.88 

6.00 

nmn 

6.22 

0.87 

Percentage 

44.0% 

44.0% 
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So  far  we  compare  MMAE  results  with  MLE  results.  Now  we  compare  MMAE 
with  all  the  other  techniques  shortly.  Tables  14  shows  the  overall  averages  of  all 
techniques  and  Table  15  shows  the  percentages  of  MMAE  statistics  equal  or  better  than 
other  techniques’  statistics  for  76  runs.  Although  there  is  not  a  best  technique  in  terms  of 
all  statistics,  AR(2)  ULS  give  more  better  results  than  others.  When  we  look  at  Table  14, 
all  statistics  of  MMAE  are  equal  or  better  than  other  techniques  statistics.  In  terms  of  ME 
MMAE  is  the  8th,  in  terms  of  MAE  MMAE  is  the  4th,  and  in  terms  of  SSE  MMAE  is  the 
1 1th  best  forecasting  technique.  The  averages  of  MMAE  do  not  seem  better  at  Table  14, 
but  if  we  look  at  the  Table  15,  MMAE  has  higher  percentages.  When  we  evaluate  both 
table  stogether,  it  is  seen  that  MMAE  is  much  better  than  all  of  the  classical  techniques. 

As  we  highlight  the  best  statistics  in  each  column,  MMAE  give  the  best, 
narrowest  3-step  and  5-step  prediction  interval  widths  with  only  3%  less  coverage  than 
average  coverage.  It  gives  the  3rd  best  l-step  prediction  interval  with  the  same  coverage 
with  the  others. 

As  a  final  word,  MMAE  forecasts  are  as  good  as  other  well-known  techniques, 
sometimes  even  better.  Especially  MMAE’s  1-step  ahead  forecast  is  better  than  all  of  the 
others.  Generally  MMAE  gives  narrowest  prediction  interval  coverages  for  all  forecasts. 
However  as  a  usual  drawback  of  these  narrow  intervals  its  coverages  are  approximately 
3%  less  than  the  other  techniques  coverages,  with  a  theoretical  coverage  probability  of 
90%. 
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Table  14.  Overall  Averages  of  Statistics  for  all  Forecasting  Techniques 
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CHAPTER  V:  CONCLUSIONS  AND  RECOMMENDATIONS 


5.1  Conclusions 

A  new  contribution  has  been  made  to  time  series  forecasting.  Our  method  is  more 
robust  and  applicable  to  any  kind  of  time  series  with  less  effort  compared  to  most  other 
methods.  We  are  saying  robust,  because  we  can  fit  MMAE  to  any  kind  of  data,  as  we  fit 
it  to  AR (p),  MA(g),  and  complex  ARMA  (3,3)  generated  data.  Even  though  we  use  only 
small  portion  of  the  possible  parameter  space,  the  MMAE  forecasts  are  as  good  as  any 
other  well-known  methods.  The  simulation  results  given  at  Tables  in  Appendix  B  and  the 
analysis  in  Chapter  IV  show  that  most  of  the  time  MMAE  statistics,  ME,  MAE,  SSE,  are 
much  better  than  most  of  the  other  methods’.  Moreover,  in  some  cases  while  some  of  the 
methods  give  poor  forecasts  when  the  data  generation  parameters  are  redundant,  MMAE 
has  never  this  problem  (see  Appendix  B  Table  98).  MMAE  can  adapt  for  changing 
parameters  and  variances. 

As  a  formal  comparison,  we  compare  MMAE  forecasts  with  ARMA  (2,2)  MLE 
(maximum  likelihood  estimation)  forecasts  and  give  the  results  in  Section  4.5.  The 
average  of  percentages  of  the  results  presented  in  Section  4.5  and  in  the  Appendix  B  is 
about  45%,  which  is  very  close  to  50%.  We  can  conclude  that  MMAE  forecasts  are  as 
good  as  MLE  forecasts,  which  is  considered  to  be  the  standard  and  hard  to  beat. 

The  best  part  of  MMAE  forecast  is:  at  the  confidence  level  of  a  =0.1,  almost 
every  time  it  gives  the  narrowest  prediction  intervals.  Although  MMAE’s  prediction 
intervals  are  narrower  than  the  others’,  its  prediction  interval  coverages  are  just  like  the 
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others’,  around  theoretical  coverage  of  90%;  the  narrower  the  prediction  interval,  the  less 
the  uncertainty  in  forecasting. 

All  these  simulation  results  and  statistics  show  that  the  MMAE  algorithm,  which 
is  taken  from  engineers,  is  worth  of  using  in  time  series  forecasting.  This  is  especially 
true  when  we  have  to  make  forecasts  frequently  on  a  daily  base  or  weekly  or  monthly 
base,  the  automated  algorithm  of  MMAE  comes  out  superior  to  others.  To  apply  other 
methods,  we  have  to  calculate  autocorrelation  and  partial  autocorrelation  functions  of 
time  series  data  and  identify  which  model  is  more  appropriate.  Later,  once  the  model  is 
identified  we  have  to  search  for  the  optimal  parameters  of  the  model,  then  we  can 
forecast.  Each  time  we  get  a  new  datum  we  just  cannot  say  that  the  old  model  is  still 
applicable.  We  have  to  go  through  whole  process  from  beginning  to  end.  As  we  talk 
about  the  advantages  of  MMAE  in  Chapter  II  and  Chapter  III,  we  say  that  we  do  not  need 
the  model  identification  and  parameter  estimation  steps.  MMAE  algorithm  automates 
those  steps  with  Bayesian  probabilistic  weighting.  We  simply  insert  the  data,  run  the 
model  and  get  the  results.  Whenever  we  get  a  new  datum,  Kalman  filter  algorithm  in 
MMAE  adjusts  to  it  and  updates  the  estimate. 

Another  benefit  of  MMAE  is  that,  while  the  other  methods  use  some  sort  of 
nonlinear  search  and  estimation  methods  (like  least  square  estimation  or  maximum 
likelihood  estimation)  to  estimate  their  unknown  parameters  [7],  MMAE  does  not  need 
any  kind  of  parameter  estimation  or  search  techniques.  It  simply  uses  all  possible 
parameters  in  the  parameter  space.  With  an  interval  of  0.1  for  parameters,  MMAE  uses 
more  than  100,000  models,  weights  them  probabilistically,  uses  the  most  likely  models 
and  gives  a  forecast,  which  is  probabilistically  weighted  mean  of  all  models’  forecasts. 
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As  a  final  verification,  we  applied  MMAE  in  time  series  forecasting  and  obtained 
good  results  in  a  Monte  Carlo  simulation.  We  recommend  MMAE  when  forecasting  is 
done  on  regular  basis.  Its  automatic  capability  allows  people  to  easily  make  forecast  as 
good  as  most  other  complex  methods,  without  knowing  anything  about  time  series 
forecasting  methods. 

5.2  Recommendations 

The  MMAE  algorithm  developed  in  this  research  is  a  viable  method  for  time 
series  analysis.  There  are  some  points  that  warrant  further  research. 

1.  We  use  a  fixed  data  mean  while  generating  data.  This  allows  us  to  forecast  only 
stationary  time  series.  However,  most  real  time  series  are  non-stationary  and  seasonal. 

We  could  build  a  model  to  estimate  the  data  mean  from  data  along  with  the  time,  for  non- 
stationary  time  series.  We  could  also  build  seasonal  Kalman  filter  to  include  in  MMAE 
bank. 

2.  We  calculate  prediction  intervals  according  to  Makridas  and  others  [15].  They  can  also 
be  calculated  according  to  Box-Jenkins  [3:142-145]. 

3.  Each  filter  in  MMAE  bank  is  started  in  a  conditional  manner,  providing  starting  prior 
values  and  shocks.  Another  approach  could  be  an  unconditional  manner,  which  requires 
backward  forecasting  with  each  filter  to  determine  the  prior  observations  before 
beginning  forward  estimation. 
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APPENDIX  A:  Data  Generation  Parameters 


Table  16.  ARMA  (2,2)  Data  Generation  Parameters 


Noise  Std.  Dev 


Case  Phi  1  Phi  2  Theta  1  Theta  2  Sample  Size  Dynamic  Measurement 


1 

0.4 

0.4 

0.2 

0.2 

300 

1 

0.1 

2 

0.4 

0.3 

-0.4 

0.3 

300 

1 

0.1 

3 

0.5 

0.3 

-0.4 

-0.5 

300 

1 

0.1 

4 

0.4 

0.5 

-0.7 

-0.5 

200 

1 

0.1 

5 

0.7 

0.2 

0.8 

-0.1 

300 

1 

0.1 

6 

0.3 

0.2 

0.6 

-0.4 

200 

1 

0.1 

7 

-0.4 

0.2 

0.5 

0.3 

200 

1 

0.1 

8 

-0.4 

0.2 

0.5 

0.3 

400 

1 

0.1 

9 

-0.2 

0.1 

-0.4 

0.5 

300 

1 

0.1 

10 

-0.1 

0.1 

-0.6 

-0.4 

800 

1 

0.1 

11 

-0.2 

0.6 

-0.5 

-0.5 

200 

1 

0.1 

12 

-0.3 

0.3 

0.8 

-0.1 

300 

1 

0.1 

13 

-0.7 

-0.6 

0.2 

0.2 

200 

1 

0.1 

14 

-0.7 

-0.6 

-0.4 

0.5 

300 

1 

0.1 

15 

-0.7 

-0.7 

-0.3 

-0.3 

100 

1 

0.1 

16 

1 

-0.6 

0.4 

0.2 

200 

1 

0.1 

17 

0.6 

-0.5 

-0.3 

0.3 

500 

1 

0.1 

18 

0.6 

-0.5 

-0.3 

0.3 

300 

1 

0.1 

19 

0.6 

-0.5 

-0.6 

-0.5 

600 

1 

0.1 

20 

0.6 

-0.5 

-0.6 

-0.5 

400 

1 

0.1 

21 

0.6 

-0.5 

-0.6 

-0.5 

200 

1 

0.1 

22 

0.3 

-0.3 

1.4 

-0.7 

200 

1 

0.1 

23 

0.4 

-0.3 

1 

-0.6 

300 

1 

0.1 

24 

0.4 

0.4 

0.2 

0.2 

300 

5 

1 

25 

0.4 

0.3 

-0.4 

0.3 

300 

5 

1 

26 

0.5 

0.3 

-0.4 

-0.5 

300 

5 

1 

27 

0.7 

0.2 

0.8 

-0.1 

300 

5 

1 

28 

-0.4 

0.2 

0.5 

0.3 

400 

5 

1 

29 

-0.2 

0.1 

-0.4 

0.5 

300 

5 

1 

30 

-0.2 

0.6 

-0.5 

-0.5 

200 

5 

1 

31 

-0.7 

-0.6 

-0.4 

0.5 

300 

5 

1 

32 

-0.7 

-0.7 

-0.3 

-0.3 

100 

5 

1 

33 

1 

-0.6 

0.4 

0.2 

200 

5 

1 

34 

0.6 

-0.5 

-0.3 

0.3 

300 

5 

1 

35 

0.6 

-0.5 

-0.6 

-0.5 

400 

5 

1 

36 

0.4 

-0.3 

1 

-0.6 

300 

5 

1 

Table  17.  AR  (1)  Data  Generation  Parameters _ 

Noise  Std.  Dev. 


Case  Phi  1  Phi  2  Theta  1  Theta  2  Sample  Size  Dynamic  Measurement 


37 

0.5 

0 

0 

0 

300 

5 

1.0 

38 

-0.5 

0 

0 

0 

200 

1 

0.1 

39 

0.5 

0 

0 

0 

200 

1 

0.1 

40 

-0.5 

0 

0 

0 

300 

5 

1.0 

41 

0.8 

0 

0 

0 

300 

5 

1.0 

42 

0.8 

0 

0 

0 

300 

1 

0.1 

43 

0.8 

0 

0 

0 

300 

5 

1.0 

44 

0.8 

0 

0 

0 

300 

1 

0.1 

45 

-0.8 

0 

0 

0 

300 

1 

0.1 

46 

-0.2 

0 

0 

0 

300 

5 

1.0 

Table  18.  AR  (2)  Data  Generation  Parameters 

Noise  Std.  Dev. 

Case  Phi  1  Phi  2  Theta  1  Theta  2  Sample  Size  Dynamic  Measurement 

47 

0.9 

-0.6 

0 

0 

300 

5 

1.0 

48 

-0.9 

-0.6 

0 

0 

300 

5 

1.0 

49 

0.5 

0.2 

0 

0 

300 

5 

1.0 

50 

-0.5 

0.2 

0 

0 

300 

5 

1.0 

51 

-0.5 

-0.5 

0 

0 

200 

5 

1.0 

52 

0.5 

-0.5 

0 

0 

200 

5 

1.0 

53 

0.5 

0.2 

0 

0 

300 

1 

0.1 

54 

-0.5 

0.2 

0 

0 

300 

1 

0.1 

55 

-0.5 

-0.5 

0 

0 

200 

1 

0.1 

56 

0.5 

-0.5 

0 

0 

200 

1 

0.1 

_ Table  19.  MA  (1)  Data  Generation  Parameters _ 

Noise  Std.  Dev. 

Case  Phi  1  Phi  2  Theta  1  Theta  2  Sample  Size  Dynamic  Measurement 


57 

0 

0 

-0.9 

0 

500 

5 

1 

58 

0 

0 

-0.9 

0 

500 

1 

0 

59 

0 

0 

0.9 

0 

500 

1 

0 

60 

0 

0 

-0.5 

0 

200 

5 

1 

61 

0 

0 

-0.5 

0 

200 

5 

1 

62 

0 

0 

0.5 

0 

200 

1 

0 

63 

0 

0 

-0.5 

0 

200 

1 

0 

57 


_ Table  20.  MA  (2)  Data  Generation  Parameters _ 

_ Noises _ 

Case  Phi  1  Phi  2  Theta  1  Theta  2  Sample  Size  Dynamic  Measurement 


64 

0 

0 

0.5 

0.4 

400 

5 

1.0 

65 

0 

0 

-0.5 

0.3 

400 

5 

1.0 

66 

0 

0 

-1.3 

-0.7 

300 

5 

1.0 

67 

0 

0 

0.5 

-0.3 

400 

5 

1.0 

68 

0 

0 

0.5 

0.4 

400 

1 

0.1 

69 

0 

0 

-0.5 

0.3 

400 

1 

0.1 

70 

0 

0 

-1.3 

-0.7 

400 

1 

0.1 

Table  21.  ARMA  (3,3)  Data  Generation  Parameters _ 

Noise  Std.  Dev. 


CasePhi  IPhi  2Phi3Theta  1  Theta  2Theta3Sample  Size  Dynamic  Measurement 


71 

0.9  -0.6 

0.5 

-0.8 

-0.5 

0.2 

200 

5 

1.0 

72 

0.9  -0.6 

0.5 

-0.8 

-0.5 

0.2 

200 

1 

0.1 

73 

0.4  0.3 

0.1 

0.7 

-0.3 

0.4 

200 

5 

1.0 

74 

0.4  0.3 

0.1 

0.7 

-0.3 

0.4 

200 

1 

0.1 

75 

0.7  -0.8 

0.3 

0.6 

0.3 

-0.8 

200 

1 

0.1 

76 

0.7  -0.8 

0.3 

0.6 

0.3 

-0.8 

200 

5 

1.0 
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APPENDIX  B:  Monte  Carlo  Simulation  Results 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  49.3  48.4 


Table  23.  Monte  Carlo  Simulation  Results  for  Case  2 
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MMAE  ARMA(2,2)  42.7  40.4 


Table  24.  Monte  Carlo  Simulation  Results  for  Case  3 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  35.1  37.5 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  55.4  57.1 


Table  26.  Monte  Carlo  Simulation  Results  for  Case  5 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  43.8  47 


Table  28.  Monte  Carlo  Simulation  Results  for  Case  7 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  53.4  53.5 


Table  29.  Monte  Carlo  Simulation  Results  for  Case  8 


^  CM  U)  CO 

o  o  o  o 


II  II  i-  CM 
03  03 

.7=  fi!  03  0 
ZZ  ZZ  JZ  ZZ 
0-  0-  1—1— 


OCDOi-OCOOO^tOOOOi-OO 

0)000)0)0)000)0)0)0)0)0)0)0)0)0) 

dddddddddddddddo 


oo)ooioiococDir)coir)ioiou)N(ON 

M-M*cococou)cocqcqcocococqooM;M; 

CMCOCMCMCMCMCviCMCMCMCMCMCMCMCMCvi 


CO  CM  N  LO  LO 


COi-NCMi-i-t-CMCMO)00 


jgqqqqqqqqqqqqqqqqq 

Qddddddddddddddddo 

E 


COCOCMLOM-lOCOCONNlOlOCDCOi-COO 

COIONCOCOCOCOCOCMM;COCOCOCOM;COCO 

CMCMCOCMCMCMCMCMCMCMCMCMCMCMCMCMCvi 


-i-CT>r^N-r^N.-»-CO^tOtOOCDOO^i- 

0)0000c00000c0c0000)c00)000)0)0)0) 

ddddddddddddddddo 


OOOCOCONNOCOM-NCOr-NOOOOO) 

COCOCM0tCOO)COONONCMCO(qSN 

^CM^CMCMCMCMCMCM^CMt^CM^^t^^ 


Ml  in  CO  CO  00  CO  CD  LO  °0 

S^cocdcM^^r^^o 


UlCDCOCON^r-OOi-N 

<OCOOOCMCMCMM‘CMr- 


CO  T-  o 

N  O  CM 


03  co 
O)  o  oo  oo 

00  O)  O)  O) 


CO  CO  1-  N  00  CM  CO 
T-  T-  CM  O  O  T- 


T-CMNi-CMi-tT-r-r-r-r-i-T-^-ri- 

ooooooqqqqqqqqqqo 

ddddddddddddddddo 


0)  a)  a)  c 

O)  D)  Q)  fc 
CO  ca  ^  CO 

5  §  ®  1 


>.  >  ■#-*  O  C/5  C/5  t/5  C/5  -  .  - 

■r  t  Q.  O  h-  1-  CNJ  CM,  "> 

OO.I  g-oDCtt<<CCCC 
2  2w  LU0C<<^^<< 


c  c 

>  >  c* 
o  o  c 


66 


PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  42  39.6 


Table  30.  Monte  Carlo  Simulation  Results  for  Case  9 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  46.1  44.2 


Table  31.  Monte  Carlo  Simulation  Results  for  Case  10 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  53.5  55.5 


Table  32.  Monte  Carlo  Simulation  Results  for  Case  1 1 
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PERCENT  OF  ESTIMATES  :R  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  29.9  26.7 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  28.7  27.9 


Table  34.  Monte  Carlo  Simulation  Results  for  Case  13 
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MAE  SSE 

MMAE  ARMA(2,2)  50.1  49.4 


Table  35.  Monte  Carlo  Simulation  Results  for  Case  14 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  37.3  37.4 


Table  36.  Monte  Carlo  Simulation  Results  for  Case  15 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  35.3  34.4 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  43.4  40.6 


Table  38.  Monte  Carlo  Simulation  Results  for  Case  17 
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Table  40.  Monte  Carlo  Simulation  Results  for  Case  19 
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Table  41.  Monte  Carlo  Simulation  Results  for  Case  20 
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Table  42.  Monte  Carlo  Simulation  Results  for  Case  21 
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Table  43.  Monte  Carlo  Simulation  Results  for  Case  22 
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Table  44.  Monte  Carlo  Simulation  Results  for  Case  23 
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Table  46.  Monte  Carlo  Simulation  Results  for  Case  25 
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Table  47.  Monte  Carlo  Simulation  Results  for  Case  26 
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Table  48.  Monte  Carlo  Simulation  Results  for  Case  27 
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Table  49.  Monte  Carlo  Simulation  Results  for  Case  28 
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Table  50.  Monte  Carlo  Simulation  Results  for  Case  29 
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Table  52.  Monte  Carlo  Simulation  Results  for  Case  31 
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Table  53.  Monte  Carlo  Simulation  Results  for  Case  32 
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Table  54.  Monte  Carlo  Simulation  Results  for  Case  33 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  42.6  41.8 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  36.1  38.4 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  45.9  45.1 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  38.9  35.7 


Table  58.  Monte  Carlo  Simulation  Results  for  Case  37 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  45.3  38.1 


Table  59.  Monte  Carlo  Simulation  Results  for  Case  38 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  54  52.4 


Table  60.  Monte  Carlo  Simulation  Results  for  Case  39 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 
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Table  61.  Monte  Carlo  Simulation  Results  for  Case  40 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 
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Table  62.  Monte  Carlo  Simulation  Results  for  Case  41 
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Table  63.  Monte  Carlo  Simulation  Results  for  Case  42 
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Table  64.  Monte  Carlo  Simulation  Results  for  Case  43 
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Table  65.  Monte  Carlo  Simulation  Results  for  Case  44 
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Table  66.  Monte  Carlo  Simulation  Results  for  Case  45 
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Table  67.  Monte  Carlo  Simulation  Results  for  Case  46 
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Table  68.  Monte  Carlo  Simulation  Results  for  Case  47 
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Table  69.  Monte  Carlo  Simulation  Results  for  Case  48 
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Table  70.  Monte  Carlo  Simulation  Results  for  Case  49 
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Table  71.  Monte  Carlo  Simulation  Results  for  Case  50 
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Table  72.  Monte  Carlo  Simulation  Results  for  Case  51 
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Table  73.  Monte  Carlo  Simulation  Results  for  Case  52 
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Table  74.  Monte  Carlo  Simulation  Results  for  Case  53 
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Table  75.  Monte  Carlo  Simulation  Results  for  Case  54 
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Table  77.  Monte  Carlo  Simulation  Results  for  Case  56 
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Table  78.  Monte  Carlo  Simulation  Results  for  Case  57 
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Table  80.  Monte  Carlo  Simulation  Results  for  Case  59 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  50.8  52 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  49  49.5 


Table  82.  Monte  Carlo  Simulation  Results  for  Case  61 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  49  49.5 


Table  83.  Monte  Carlo  Simulation  Results  for  Case  62 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  49  51.5 


Table  84.  Monte  Carlo  Simulation  Results  for  Case  63 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  55.3  55.9 


Table  85.  Monte  Carlo  Simulation  Results  for  Case  64 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  37.7  38.5 


Table  86.  Monte  Carlo  Simulation  Results  for  Case  65 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  42.5  44.1 


Table  87.  Monte  Carlo  Simulation  Results  for  Case  66 
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MAE  SSE 

MMAE  ARMA(2,2)  41.3  42.8 


Table  88.  Monte  Carlo  Simulation  Results  for  Case  67 
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MAE  SSE 

MMAE  ARMA(2,2)  35.6  37.6 


Table  89.  Monte  Carlo  Simulation  Results  for  Case  68 
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MMAE  ARMA(2,2)  48.8  50 


Table  90.  Monte  Carlo  Simulation  Results  for  Case  69 
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MAE  SSE 

MMAE  ARMA(2,2)  51.2  54 

MMAE  ARMA(2,2)  85.6  85.6 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  54.8  51 .2 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  47.6  49.2 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  53.6  52.8 


Table  94.  Monte  Carlo  Simulation  Results  for  Case  73 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  49.2  44.8 


Table  95.  Monte  Carlo  Simulation  Results  for  Case  74 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  52.4  52.8 


Table  96.  Monte  Carlo  Simulation  Results  for  Case  75 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  24.4  24.4 


Table  97.  Monte  Carlo  Simulation  Results  for  Case  76 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  25.6  27.2 


Table  98.  An  Example  of  Monte  Carlo  Simulation  Results  Showing  Robustness  of  MMAE 
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PERCENT  OF  ESTIMATES  BETTER  THAN  MLE  OUT  OF  1000  REPLICATIONS 

MAE  SSE 

MMAE  ARMA(2,2)  85.6  85.6 
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